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Abstract. Amongst the bound states produced by the strong interaction, radially excited 
meson and nucleon states offer an important phenomenological window into the long- 
range behavior of the coupling constant in Quantum Chromodynamics. We here report on 
some technical details related to the computation of the bound state’s eigenvalue spectrum 
in the framework of Bethe-Salpeter and Eaddeev equations. 


1 Introduction: excited states as an eigenvalue problem 

A great deal of research activity in hadron physics is concerned with hadron structure and revolves 
around two fundamental questions: what constituents are the hadrons made of and how does Quantum 
Chromodynamics (QCD), the strong interaction component of the Standard Model, produce them? 
These are simple questions which, however, may not entail simple answers. To understand the mea¬ 
surable content of QCD, spectroscopy is a valuable and time-honored tool — suffice it to mention 
the inestimable progress made in the computation of atomic or molecular spectra and subsequent 
comparison with experiments that lead to a deeper understanding of Quantum Electrodynamics. 

The same is true for QCD: if confinement is related to the analytic properties of QCD’s Schwinger 
functions, then light-quark confinement ought to be understood by mapping out the infrared behavior 
of the theory’s universal fi function. Obviously, this cannot be possibly achieved in perturbation 
theory. A nonperturbative continuum approach to QCD is provided by Dyson-Schwinger equations 
(DSE) which relate the theory’s fS function to experimental observables m. Therefore, comparison 
between DSE predictions embedded in bound-state calculations and observations of the hadron mass 
spectrum as well as of elastic and transition form factors can be used to study the long-range behavior 
of QCD’s interaction. As the properties of excited hadron states are considerably more sensitive to the 
long-range behavior of the strong interaction than those of ground states nan, excited mesons and 
nucleons lO are an important source of information and complement our understanding of the strong 
interaction in light mesons 00 and in heavy-light mesons with disparate energy scales nsHsi- 

Of course, the properties of these hadrons can only be understood, and the functional behavior of 
the p function be inferred therefrom, by studying the quark’s DSE in conjunction with quark-antiquark 
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or three-quark bound-state equations, the Bethe-Salpeter and Faddeev equations, respectively. Both 
are treated as eigenvalue problems which we here exemplify with the relativistic bound-state equation 
of pseudoscalar j’’ - 0“ mesons. The homogeneous Bethe-Salpeter equation for this qq bound state 
with relative momentum p and total momentum P can be generally written as, 

ro-(p,P) = J ^{pXP)[S(k + q^P)Yf,-{k,P)Sik-q^P)\ , (1) 


where Sik ± q±P) are dressed quark propagators with q+ + - 1 (NB: in a Poincare invariant 

calculation numerical results are independent of the momentum partition parameter q±). For the sake 
of simplicity, we omit flavor and Dirac indices, as the following discussion is independent of them. In 
the same spirit, we restrict ourselves to the rainbow-ladder truncation of the interaction kernel. 


'K{p,k,P) = - 


Zlgiq^) 


A‘‘ 

-jJfiT^yiq) —Jy , 


( 2 ) 


with the transverse projection operator T^y{q) g^y - q^qyjq^, q - p - k, Z 2 is the wave-function 
renormalization constant and A" are the SU(3) color matrices in the fundamental representation. Var¬ 
ious model ansatze iT4l4T7]l have been proposed for the effective interaction, §{q^), which emulates 
the combined effect of the gluon and quark-gluon vertex dressing functions and most recent efforts 
extend these models to include important transverse components of this vertex beyond the leading 
truncation US). Most importantly, Eq. Q satisfies the axialvector Ward-Takahashi identity 1(191 and 
therefore ensures a massless pion in the chiral limit. As we shall see below, Eqs. Q and Q define 
an eigenvalue problem with physical solutions at the mass-shell points, P^ - -m, where m is the 
bound-state mass. 

The Bethe-Salpeter equation’s Poincare-invariant solutions can be cast in the form. 


To- {p,P) ^ ys[ilDEQ-{p,P) + rFQ-{p,P) + p{p-P)Go-{p,P) + (T^yp,PyH^-{p,P)] , (3) 


Note that this Euclidean-metric basis, JT’{p,P) - ysiHD, f, Pip ■ P), cr^yp^Py], is nonorthogonal 
with respect to the Dirac trace. The functions T^Xp, P) - [Eo-ip? P)^ Fo-(p, P), Go-(p, P), Elo-(p, P)} 
are Lorentz-invariant scalar amplitudes and are extracted from the Bethe-Salpeter amplitude ([^ with 
appropriate projectors, P^ip, P)|^ 

i P) Ti-d [xf(p, P)XPip, P)] = day ; P"(p, P) = P“^(P, P) Xf{p, P) , (4) 

where 0,(5 - 1, ...,4 and the projection is given by: 

Po“ ip, P) = 3 Tro [pXp, P) To- ip, P)] ■ (5) 

Using the Bethe-Salpeter equation, this leads to the eigenvalue problem. 


A(P^) (p, P) = k, P) (k, P ), 

where 'K°^^{p, k, P) stems from the projection of Eq. ([^ using Eqs. (j^ and (j^: 

'K^PipXP) = -^^^T^y(q)\xcoo[P“ip,P)y^,A^Siky)X{f^ik,P)SikS)yyA‘‘] . 


( 6 ) 


(7) 


^ see, e.g.. Ref. ESI for the derivation of the change of basis coefficients P). 
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In Eq. A(P^) is a scalar function and the eigenvalue equation has a solution for every value of 
Typically, iterative eigenvalue algorithms are employed as only under simplifying assumptions an 
inversion of the Bethe-Salpeter problem is possible lEIHil- To elucidate the iterative procedure, we 
simplify Eq. (j^ (see also Ref. Il24l l. so that; 

A(F^) lO) = Id)). (8) 

The kernel, 'K(P^), has a complete set of real eigenvectors with eigenvalues Ai{P^) which are 
ordered as AaiP^) > A\(P^) > A 2 (P^) > .... > AAP^). Thus, any solution can be written as a linear 
superposition, 

OO 

= (9) 

i'=i 

where a, are real constants and the vector must be normalized. To begin the iterative process, one may 
“guess” as solution, such as in Eq. (j^, for a given value of P^ and n successive actions of the kernel 
lead to, 

OO 

1=1 

Since Aq > A/, the coefficients of 10,) converge to zero for sufficiently large values of n and therefore 
the amplitude | converges to the ground state amplitude | (po): 

|0„>"=“/lSflol0o) - /lo-7C"-'(p2)|<i)). (11) 

This is the most basic method of computing the largest eigenvalue Ao{P^) and its associated eigen¬ 
vector to any required accuracy and is referred to as power or von Mises iteration. The trajectory of 
this eigenvalue function for a range of P^ values yields the ground-state meson mass, that is one finds 
Ao(P^) - 1 when P^ - -tn^. 

2 Gram-Schmidt orthogonalization 

A common procedure to extract the wave function of excited states from the spectrum of the inter¬ 
action kernel is based on the Gram-Schmidt orthogonalization process. As mentioned in Section [T] 
the homogeneous Bethe-Salpeter equation is an eigenvalue problem, where the largest eigenvalue the 
interaction kernel produces corresponds to the ground state. However, the eigenvalue spectrum is not 
limited to the ground state and excited states with smaller eigenvalues can be determined with the 
same iterative methods discussed above. As sketched, e.g. in Ref. Il24l . one chooses a mass that is 
larger than the ground state mass, P^ - -m^ < -m^, and applies the iterative procedure outlined in 
Section[T] 

• find the largest eigenvalue, To, and the associated eigenvector for and /Io(ot^) > /lo(OTg). 

This is the unphysical “ground state” at the mass scale P^ - -m^. 

• make again a guess for the Bethe-Salpeter amplitude now projecting out the eigenvector that pertains 
to the largest eigenvalue /lo(m^) in a first iteration. 

• use the eigenvector obtained in the previous step as input for the Bethe-Salpeter amplitude and 
project out again the eigenvector associated with Ao(irP') in a second iteration; the resulting eigen¬ 
vector must be projected as before in a third iteration and so on. 

• the iteration converges after the nth projection which yields an eigenvalue, Ai(m^), with an eigen¬ 
vector orthogonal to the one associated with the “ground state” eigenvalue Aoinp-) > Ai(m^). 
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Figure 1. First two eigenvalue trajectories of the quark-diquark kernel used in Refs. 12511261 as a function of 
V-7^. Note that the eigenvector solutions are here differently normalized so that the eigenvalues are d, = 1 lg\ 
instead of d, = 1 for see Ref. (25) for details. The intersections of the solid horizontal line with the 

trajectories locates the mass-pole position of the nucleon and its first excited state identified with the Roper. 


• by varying one obtains the mass evolution of the second largest eigenvalue A\(P^), which is 
exemplified in Figure The solution for Ai(P^) = 1 yields the mass of the first excited state, 
P^ - -OTj where mi > mo, at which the eigenvector is the state’s Bethe-Salpeter amplitude. 

The normalized projection within this Gram-Schmidt procedure is effected by, 

<00 I <E)) 

I®) - I®) - 77 ^ I0o), (12) 

Wo 100 ) 

where lO) is the initial guess for the Bethe-Salpeter amplitude in Eq. ([^ and | 0o) is the vector of the 
ground-state Bethe-Salpeter amplitude. In projecting out this ground-state, one must define a norm in 
Euclidean space via an inner product, as is evident from Eq. ( [T^ . Within Quantum Eield Theory in 
rainbow-ladder approximation, this product is defined as: 

<'I'|(l>>:=trcDj ^ik,-P)S{k + jj^P)<S>(k,p)Sik-Ti^P). (13) 


Orthogonality is thus defined by the condition, <'I'|<I)) = 0, which expresses the vanishing of the 
overlap amplitude at P^. Note that Eq. (13 1 is not valid for computations beyond the leading approxi¬ 
mation, i.e. the rainbow-ladder truncation. 


3 Krylov subspace and Arnold! iteration 

After momentum discretization and Chebyshev expansion of Eq. (|^, the numerical kernel of the 
Bethe-Salpeter equation (or the Eaddeev equation for the three-body problem) is a non-symmetric 
matrix of large dimensions and with some of its eigenvalues close to zero. In general, non-symmetric 
matrices have eigenvalues and eigenvectors which are very sensitive to small changes in the matrix 
elements due to the lack of symmetries on which traditional methods rely to ensure numerical stabil¬ 
ity 127]. This is particularly true for radial excitations, where angular dependence encoded in higher 
Chebyshev moments contributes in a nontrivial manner. Therefore, we make use of the numerical 
ARPACK library ll28ll which is designed to compute a few eigenvalues and corresponding eigenvectors 
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of a general nxn matrix by means of the implicitly restarted Arnold! method (IRAM). The strength 
of the Arnold! method Il29l lies in the application of the stabilized orthogonalization algorithm in the 
Krylov space of a given matrix K to find its eigenvectors, as will be explained shortly. The most 
important feature of this algorithm is its ability to decouple the calculation of the eigenvectors corre¬ 
sponding to the eigenvalues with the largest absolute value from the eigenvectors whose eigenvalues 
have an absolute value close to zero. This is highly desirable as the calculation of the latter are plagued 
by numerical instabilities. 

In algebra, the order-r Krylov subspace is generated by an operator whose representation is given 
by an n X n matrix K - Kij, and a vector <1) of dimension n. It corresponds to the linear subspace 
spanned by actions of K on O: 

Sr := '<!>} . (14) 

The power iteration of K yields the sequence in Eq. and converges to the eigenvector associ¬ 
ated with the largest eigenvalue do. However, hereby one only makes very limited use of the stored 
information, as only the final vector, ’<E), is kept. On the other hand, the basis for the Krylov 
subspace is derived from the Cayley-Hamilton theorem which implies that the inverse of a matrix 
can be expressed as a linear combination of its powers. Thus, Krylov subspaces play an important 
role in contemporary iterative methods to obtain one or few eigenvalues of large sparse matrices or 
to solve large systems of linear equations. Instead of heavy matrix operations, these methods rely on 
successive multiplications of vectors by the matrix, thus forming a Krylov subspace, and then employ 
the resulting vectors. 

The vectors of the Krylov space are initially not ortogonal and usually become almost linearly 
dependent due to the properties of the matrix power iteration. Nonetheless, an orthogonal matrix 
can be constructed from the basis vectors by means of the Gram-Schmidt process described in Sec¬ 
tion 1^ This method proves to be unstable but its shortcoming can be overcome with the Arnold! 
iteration ||2^ which uses the stabilized Gram-Schmidt process and can be applied to general, possibly 
non-Hermitian matrices. The Arnold! method generalizes the Gram-Schmidt process by computing 
the eigenvalues of the orthogonal projection of K onto the Krylov subspace, where the projection is 
represented by the upper Hessenberg matrix Hr ll28l . For Hermitian (symmetric) matrices, the Arnold! 
iteration is analogous to the Lanczos iteration. 


4 A note on orthogonality 

In the context of Bethe-Salpeter and Faddeev equations, orthogonality is defined by Eq. ( [T3] l. In 
vectorial form, using the shorthand x - and y - p^, Eq. (j^ can be written as. 


4(l/) FR(y) ^ Jdx 7C(r, y) ■ Fr(x) , 
which in a numerical treatment is evaluated as the sum, 

A(xj)FR{xj) = ^ Wi'K{xj,Xi) ■ FR(xi) , 
i 

and similarly. 


A{xj) Fiixj) = ^ Wi Fiixi) ■ '7C(xy, x,), 


(15) 


(16) 


(17) 


where x, are the nodes of a given quadrature with weights w, . Since the integral in Eq. 0 is over a 
charge-conjugate Bethe-Salpeter amplitude, -P) := -P)C^, we introduce in Eqs. (16 1 
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and([^ “left” and “right” eigenvectors, Fi and Fr, respectively. Multiplying Eq. ( |T6 ] i from the left 
by Tjj we obtain, 

A ^ Wj F'^ixj) ■ Ffi(xj) = ^ Wi Wj F'^ixj) ■ 'K{xj, xd ■ FR{xi) = ^ Wi F'^ixj) ■ FR(xi) , (18) 

j i.j ' 

where F'^{xj) is the eigenvector associated with the eigenvalue A'. For A + A! this relation implies, 

fdx F'^(x) ■ Fr(x) ^ 0 , (19) 


The left eigenvector is given by the Bethe-Salpeter wave function, S(k- t]-P) -P)S (k + t]+P), 
and the right eigenvector is the Bethe-Salpeter amplitude ll^ . We thus deduce from Eq. (|T^, 

trcD^5(/t-77_P)O^,(;t,-E)5(/t-H77+P)(l)^(;t,E) = 0, for/^t/l, (20) 


and because the trace is cyclic this orthogonality condition is equivalent to that described in the para¬ 
graph below Eq. ([13 1 . Moreover, we verify that the Bethe-Salpeter equation spectrum of Eq. is 
equal to that in Eq. (|l7|i. We stress that Eq. (19 1 is generally valid and no assumption was made about 
the kernel structure or about the eigenvectors. 


5 Examples in hadron physics 

The meson and nucleon resonance structure has been the object of a long history of studies and we 
here limit ourselves to the approaches based on the combined Dyson-Schwinger and Bethe-Salpeter 
(or Faddeev) equations, in particular their application to excited mesons and nucleons to compute 
their masses, weak decay constants and/or electromagnetic form factors making use of the techniques 
described in Sections [T}@] 

Seminal studies on the ground state spectrum of light pseudoscalar mesons established that the 
7r(1300) can be described as the first radially excited state of the Goldstone boson and furthermore 
that the decay constant of excited states, P„, vanishes identically in the chiral limit OTll . 

/“;°(A<) = 0,n>l, (21) 

where m is the renormalization-group invariant current-quark mass. Electromagnetic properties of 
ground and excited state pseudoscalar mesons, making use of the the Gram-Schmidt process illus¬ 
trated in Section [^ were studied in Ref. El. A hrst approach to computing the eigenvalues of 
unflavored light and heavy mesons and their corresponding Bethe-Salpeter amplitudes with the im¬ 
plicitly restarted Arnold! factorization, as implemented in the ARPACK ll28l library, is expounded 
in Refs, ll^ 1^ and has successively been applied to the light llTSl and heavy f35l quarkonia 
spectrum for pseudoscalar, vector and tensor states. A recent detailed analysis of the Maris-Tandy 
interaction m parameter space for ground and radially excited states is presented in Ref. lEl, 
where the authors conclude that the preferred inverse effective range of the Maris-Tandy interaction in 
rainbow-ladder truncation, a>, is shorter for heavy quarkonia than for the light quarkonium spectrum; 
coqq - 0.7 GeV vs. cjgg - 0.5 GeV. 

The same library has then been employed to study radial excitations of flavor-singlet and flavored 
pseudoscalar mesons within the framework of the rainbow-ladder truncation ||4l[36l|371. A summary 
of the theoretical values for the pseudoscalar’s ground- and excited-state masses and weak decay con¬ 
stants for flavor-singlet and nonsinglet = 0“ mesons is reproduced in Table for two parameter 
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Table 1. Masses and decay constants for flavor singlet and nonsinglet = 0 mesons; see 


Section]^ 



Model 1 [GeV] 

Model 2 [GeV] 

Reference 

m„ 

0.138 

0.153 

0.139 113^ 

fn 

0.139 

0.189 

0.13041381 


0.990 

1.414 

1.30 ±0.10 mi 

Ac 1300) 

-1.1 X 10-^ 

-8.3 X lO-"* 


niK 

0.493 

0.541 

0.493 11381 

h 

0.164 

0.214 

0.156 El 

'«ir(1460) 

1.158 

1.580 

1.460 El 

/rr(i460) 

-0.018 

-0.017 


m-ss 

1.287 

1.702 


f-ss 

-0.0214 

-0.0216 


i^nAiS) 

3.065 

3.210 

2.984 11381 

fndis) 

0.389 

0.464 

0.395 El 

’^ncQ-S) 

3.402 

3.784 

3.639 El 

ftU2S) 

0.089 

0.105 



sets of the interaction introduced by Qin et al. im. Namely, Model 1 and Model 2 correspond to the 
interaction parameters, a> - 0.4 GeV, toD = (0.8 GeV)^ and a> — 0.6 GeV, toD = (1.1 GeV)^; see 
Ref. El for details and Ref. El for a discussion of the effective interaction. Note that we explored 
a range of parameter combinations for loD + const., yet could not find a unique set which produces 
theoretical observables that compare well with all data in Table We therefore observe that we are 
unable to consistently and simultaneously describe both ground and excited state observables with 
the given interaction ||2l in the rainbow-ladder truncation — whereas one parameter set successfully 
reproduces experimental ground state masses and decay constants but not the experimental data on 
excited states, the second parameter set only provides a reasonable description of the excited pseu¬ 
doscalar mass spectrum. 




Figure 2. Left panel: Lowest Chebyshev moment, “Fp,(p^), associated with the leading Dirac structure Ep^{p^) 
of the pseudoscalar’s Bethe-Salpeter amplitude for the first radial excitations Pi = ;r(1300), (ii)„=i and 
Right panel: First three normalized Chebyshev moments, ‘S i(p^), i = 0,1,2, of the leading S'-wave component 
in the nucleon’s first excited-state Faddeev amplitude. 
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As it becomes clear from the left graph of Fig. for > \ GeV^ the amplitude’s lowest 
Chebyshev projection, '^Ep^ (p^) associated with the leading covariant ijs of the pseudoscalar’s Bethe- 
Salpeter amplitude, becomes negative definite in case of the first radial excitations, 7r(1300) and 
(Si)„=i, whereas it remains positive for the ground states. It is remarkable that the behavior of the 
Bethe-Salpeter amplitudes parallels the pattern of wave functions in quantum mechanics, namely, the 
number of zeros can be associated with a principal quantum number n. As just mentioned, the ground 
state amplitude has no zeros and can thus be associated with n - 0. The amplitude of the next highest 
mass mesons possesses one zero and one assigns the quantum number n = 1 and so on. An analo¬ 
gous behavior is found for the nucleon’s first excited state Faddeev amplitude for which we plot the 
first three Chebyshev moments of its leading S -wave component. Indeed, very recently a range of 
properties of the proton’s radial excitation was predicted which strongly suggests that the nucleon’s 
first radial excitation is the Roper resonance ESI . In particular, in Fig. 3 of Ref. ES\ one can ap¬ 
preciate that the overlap amplitude between the nucleon and its first excited state described by the 
nucleon-Roper transition form factor F*{Q^) vanishes at 2^ = 0 and therefore satisfies orthogonality. 
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